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ABSTRACT 


Any change in neutronic properties in a reaetor oper t- 
ing at steady state will result in a change in the equilib- 
rium neutron flux and hence, the power of the reactor. A 
main cause for a change in neutronic properties is the high 
temperature attained in a reactor, which produces a feedback 
in the reactor operation. The response of the reactor to 
a particular feedback is analyzed by using Liapunov's 
cond Method to specify stability regimes. Both, the 
point-kinetics model and the distributed parameters system 
are analyzed. 

Data from a typical thermal and fast reactors is used 


to specifically determine the stability domains. 
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I. LIAPUNOV'S SECOND METHOD 


As HLS TORY 

In 1892 the Russian mabhenatician A. Μ. Liapunoew pab- 
lished a long paper dealing with the problem of stability 
ofemotion- (1).- This paper was translated into Fremch in 
1907 and reprinted in America in 1949, Liapunov's Theory 
received by then only little attention due to the diffi- 
culty in understanding the advanced mathematical theorems, 
Ba the abstracteway itewasepresentedeand to its lack of 
practical application and for a long time it was nearly 
forgotten. About 35 years ago, Soviet mathematicians re- 
Sumed the investigation of Liapunov's Theory and its excel- 
iSt application in several technical fields, mainly in 
control engineering, was ne ded. Significant work in this 
area was published by Malkin (2), Letov (3), Lur'e(4) and 
Chetaev (5). The excellent paper by Massera (6) and the 
translation to English of most of the Russian works stirred 
up considerable curiosity in this country. Bellman (7) in 
E55 published an excellent work concerning stability, but 
Huc section on Liapunov's method is difficult to comprehend. 
According to this author, Hahn (8), Krasovskii (9) and 
LaSalle and Lefschetz (10) have the best treatises on the 
subject, the consolidation of the concepts of Liapunov's 
Theory and the clear presentation of it, make these three 


books the main references for workers in the field. 
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Recently, Zubov. (11) found the best extension of the 
Diapunow'su Theorysssto incluüudesthe analysis οἳ "σε ο 
ferential equations, which is by now the main topic for 
research. Yet, a good amount of work remains to be done. 
A the previous works, including Liapunov's original 
theory, were devoted to the analysis of stability of 
ordinary differential equations. 

Many authors have applied Zubov's extension to works 
concerning vibrations, reactor physics, hydrodynamics, 
magnetohydrodynamics and control processes. The list of 
references will provide the interested reader with the 


mon works in this fields 


Dr THEORY 

ine nane second method (or ‘direct Method is of 
maestorical origin. Liapunov s Theory is divided in two 
categories. The "first method" which comprises all pro- 
cedures in which the explicit form of the solutions is 
used, specially when represented by infinite series and the 
en eSstigation of stability “in the small" (local stability) 
Studies the singular points of a nonlinear differential 
equation by using the appropriate linearized version of the 
differential equation near the singular point. The "second 
method" attempts to make stability statements by using (in 
addition to the differential equations) suitable functions, 
called Liapunov functions, which are defined in the motion 


space. This method which deserves special study due to 
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its inherently advantages over most Sof the conventional 
methods, investigates the stability “dm the large” iom 
Stability). 

Liapunov's Second Method is in essence a more general 
dn ession of the Hamiltonian (total energy), and it is 
based in the statement that a physical system loses poten- 
tial energy in a neighborhood of a point of stable equilib- 
mum. Ilt is said that it is more general than the total 
Enerov method, because unlike the energy of a system, the 
Ex Dpunov function, denoted V, is not unique. This is the 
main reason why the "second method" is a tool in the analy- 
sis of stability of dynamical systems, which must be used 
with considerable skill. One of the main features of this 
Method is its appeal to geometric intuition; V(x,t) can be 
seen as a measure of the "distance" of the state (x,t) from 
the origin, in the state space and the variation of this 
Edeestance™ (norm of the differential equation) as t varies 
will provide definite bounds of stability regions for the 
Prescribed system under consideration. 

Normally stability with respect to one norm does not 
imply stability with respect to another norm. This diffi- 
culty does not appear in finite-dimensional systems since 
all norms defined on a finite-dimensional vector space are: 
equivalent. 

l. Fundamental Concepts and Definitions 

Stability Τα a property of certain systems of dif- 


ferential equations and is basically concerned with the 
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question of whether or not the dynamic system will return 
to a particular state after it has been disturbed in some 
way. 

Differential equatzons in-therr most gener ia 


can be expressed as 


DUE) - Plu (, c) t] 


where U(x,t) is a multidimensional function of space and 
time. It may well happen that F depends upon U(t) alone 
and not upon time explicitly. Then the previous equation 
assumes the simpler form 
Da) = 
A system of this nature is known as autonomous. Since it 
is not intended in this work to expose the reader with the 
Mathematical aspect of Liapunov's Theory, in general, the 
presentation of the theory will be made without proofs. 
Let us denote Q(R) the spherical region ||x|| <: R -and by 
A(R) the sphere boundary ||x|| = R. The matter of concern 
is the stability of the origin. Initiating the motion at 
a point Χο» the origin is said to be: 
(a) Stable whenever the path remains in the spherical 
region N(R); that is,the path never reaches the 
boundary sphere A(R), Figure 1. 
(b) Asymptotically stable whenever it is stable and 
in addition the path tends to the origin as time 


increases indefinitely. 
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Figure 1. Domains of stability. 
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(c) Unstable whenever the path reaches the boundary 
sphere A(R). 

For autonomous systems: 

A Liapūünov- function: V(U), is defined as air eie 
Me definite scalar function with the following properties: 

(a) VU) *. continuous with its- first partial deriva- 
tives in a certain pen =repionsabeut the oxigdn 

(b) V(o) » 0 only 

(c) Outside the origin and always inside the open 
region, V is non-negative and vanishes only at the 
origin. Ὅμο ορ. πο αρ an isolated minimum ol VW; 
and in addition 

V Sp Ue trees open mee ion 

If V = 0, the stability is neutral 

If V « 0, the stability is asymptotic. 

For nonautonomous systems: 

Let us define V, (U) as previously. A Liapunov 
function, V(U,t) is defined as a positive definite scalar 
function with the following properties: 

Ca) V(U,t) is defined in the open region for all 

Ea 


Cb) ο) το στο ο ο 


|^ 


( c) V (CU) VCU,t) for all x in the open region and 
ait t 270, (and jing addi tion 


τ the open region 


το 





where 





Summarizing hese conditions, he three main 
Staps to «eet if «aeccadar fupetían is a Liapunov func- 


morare: 
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28 Methods to Construct a Liapunov Function 


Several authors have attempted to present guide- 
lines for generating Liapunov Functions. Unfortunately, 
most of these methods are extremely restrictive, limita- 
tions being imposed primarily by the number of nonlinear 
terms and the order of the system.  Ingwerson (12) pre- 
sents a Table to generate Liapunov Functions for linear 
autonomous differential equations up to the fourth order. 
Barbashin (13) and Lur'e and Rozenvasser (14) tried to 
establish some rules for construction of Liapunov Functions. 
Schultz (15) has one of the best methods available now; it 
is called the Variable Gradient Method, being the least 
restrictive for the case of autonomous systems. Its only 
restriction is that all nonlinearities must be single- 
valued. 

In general, literature late in 1961 and 1962 became 
saturated with methods for generating Liapunov Functions, 
most of them with the restrictions previouslyəmentioned. 
Furthermore, such methods, when applicable, were directed 
to autonomous systems only. 

Although development of Liapunov stability theory 
and applications to ordinary differential equations has 
Progressed rapidly, its application to partial differential 
equations has remained limited. The importance of partial 
differential equations in the fields of reactor physics, 
hydrodynamics, control processes, etc. has motivated inves- 
tigations of possible ways to extend Liapunov stability 


+ 


νο: tO partlal differential equations. 
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τὸ τας study Functional Analysis plays an important 
role, because the stability of the equilibrium solution is 
defined in terms of the norm induced» by the inner product of 
the Hilbert space on which the solutions of the system are 
defined. Thus, with proper choice of inner product or norm, 
the square of the norm becomes the Liapunov functional which 
estebtirshes ecesympieticestebüli«y.  lno£ect what is being 
done is the construction of a scalar function of the distance 
between the solution and the equilibrium point of interest. 
If this distance function, evaluated along the solution 
paths, obeys certain inequalities, then M ET concern- 
ing the stability of the equilibrium point can be made. 
E Dan (16) defies stability in terms of two metrics, 
rather than one, to be more restrictive on the initial 
states. The choice of the initial dis Space and the 
Metric is crucial in the formulation of stability problems 
in the framework of Liapunov's Second Method. For cases in 
which only the behavior of some of the state variables is 
of importance or some function of the state variables is of 
interest to the analyst, then, for these cases, it is mean- 
mmetul to define Liapunov stability with respect to two 
metrics. For instance, one may be only eee im tehe 
behavior of the maximum deflections in an elastic system, 
regardless of the velocity and acceleration involved into 
attained it. 


According to Kastenberg & Ziskind(17) extreme care 


has to be taken when interpreting the obtained results, 


το 





Because severe peaking in some of rhe state vazıablesgarzene 
system can cause the L, norm oí the system trajectory to 
move out of the domain. Wang (18) and Parks (19) have de- 
voted special attention to the study of panel flutter which 
represents a linear distributed parameter system. It is 
Concluded that a good comprehension of Hilbert and Banach 
Spaces will provide the analyst {with an excellent tool for 
generating Liapunov functionals. 

3. Comparison with Other Methods 

The advantages of Liapunov's Second Method over the 
conventional methods used for stability ee can be 
summarized as follows: 

(a) It employs the systenwegqgwatjons wirt howeswesorting 
to approximations. Normally a distributed system is approxi- 
mated by a lumped parameter model having finite or infinite 
number of degrees of freedom. This method is not satisfac- 
tory because qwite oftem it exhibitsechamaetermsties which 
do not agree with the physical nature of the problen, 
Liapunov's method is in general a more reliable method. 

(b) There is no theoretical limit on ene number of non- 
linearities or on the order of the differential equation to 
which Liapunov's Second Method can be applied. 

(c) The laborious work implied in finding the system 
solution is avoided. Normally the system equations are in- 
tegrated numerically for some given perturbations in the 


inıtaalesconditions, This method is time consuming and 
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sometimes does not give an accurate presentation of the 
system behavior. 

(d) The relationshipmbetween the system staba | uw and 
the system parameters is directly extracted from the analy- 
SiS of Liapunov's methed conditions for stabmuxl ise. 

(e) There are no mathematical restrictions with res- 
Bece Co unique. 

The method presents also some deficiencies: among 
them the main ones are: 

(a) The construction of the required Liapunov Function 
for the system under consideration is without doubt the 
most difficult part of the task due to the fact that there 
mano guideline to construct it and success depends mainly 
upon the ingenuity and experience of the analyst. 

(b) The interpretation of the obtained results has to 
be done carefully. 

ο intecray inequalvties play an important role in 
deriving conditions for stability. As a consequence, the 
regions of stability may be somewhat looser than those ob- 
tained by other methods. 

Inweoncedusronjg"in spite of ehe difficulties cited 
above, the advantages of the Liapunov's Second Method makes 


it an excellent method for studying stability. 


C. APPLICATION TO NUCLEAR REACTOR CONTROL 
Liapunov's Second Method has been applied with success 


by Hsu (20) who considers the linear and non-linear cases, 
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analyzing them by means of the Spectral Analysis andas 
Lov s Method, and comparing the results obtained by omb 
method. Hsu shows that the Spectral Analysis provides a 
Beetle more information than Liapunov's Method, but this 
latter eliminates the calculation of the eigenvalues of the 
system, and for the nonlinear system analysis bounds of 
stability arc presented. 

Kastenberg (21) presents a clear comparison of Liapunov's 
Second Method with the Semigroup Method, and the Comparison 
Function and Marian Principle Techniques. In the Liapunov 
and Semigroup methods, one must obtain an a priori bound on 
the system nonlinearity and then proceed. When employing 
the Maximum principle or the comparison function technique, 
an a priori bound on the system nonlinearity is not always 
απ ττθδά. In contrast, one must give a bound on the initial 
condition. For cases which are stable in a global sense, 
ος results of the various methods coincide. 

The application of Liapunov's Second Method to nuclear 
reactor control seems to be excellent, mainly for the spa- 
tially-dependent reactor system, due to the fact that the 
method allows the inclusion of any number of nonlinearities. 
Pris will permit the study of spatial effects of temperature, 
Ate main feedback effect in a reactor, control rod motion 
and several other processes within a reactor. Also, the 
inclusion in the analysis of the different reactivity coef- 
ficients, as doppler resonance broadening and structural 


expansion is allowed. 


zu 





Às pointed before, the anterpretation cf thermo nee 
must be done very carefully due to the fact thar ar rapuney 


Eunctoon is not unique. 
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Ti. MATHEMATICAL MODEL OF- THE REACTOR oX TN 


Both Thermal and Fast reactors are described essential- 
ly by the same basic dynamic principles, regardless of 
material and geometry considerations. 

For a given reactor size, the reactor reactivity depends 
on the neutron cross-sections and on the relative amounts 
and densities of different materials. All of these being 
affected by the temperature, the reactivity will be strong- 
ly coupled to the power of the reactor and the reactor 
governing equations become nonlinear. 

The Thermal Reactor stability is analyzed by means of 
the lumped-parameter (point-kinetics) model, in which the 
partial differential equations are reduced to ordinary dif- 
ferential equations via spatial discretization. This is 
justified only when infinitesimal small perturbations are 
introduced, as in the case when k is very near to unity 
Mery small departures from "critical"). Thermal reactors 
are generally smaller in size than fast-breeder reactors. 
The Fast-Breeder Reactor stability is analyzed using the 
governing partial differential equations without resorting 
to any type of approximations. Due to possible stronger 
space effects in fast reactors than in thermal reactors, 
the space and time of the state variables of the system 
must be maintained during the analysis. The reactivity in- 
ο in the system, Figure 2, can be: either ρου ο σαν 


üegative. Both cases will be considered in this work. 
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owe POINT-KINET ICS MODEL 
l. Diffusion Equation 
According to Meghreblian and Holmes (22). ΕΠΕ ΠΡ 


muon diffusion equation is written in the following way: 


D V^$(r,t) - E = Ll ΑΣ - 


eV ρε σα ο νο ο Ὁ (1) 


di 
Mepressing the equation for a slab reactor and using only 


eme-delayed neutron group: 


p3*6Gut) a a _ 
V ot 
Ox 
- (1-8)v E, pg o(x,t) - pg A C(x,t) (2) 


It is looked for a solution which is separable in time and 
space: 

O(x,t) = g(t) F(x) | αν ον 
Ct) εἴα) (b) 


ος τε) 


Applying relations (3) to Equation (2) yields 


, ] 
moe o F(x) _ 2v. TO. a 
F(x) ,,2 u : D 
- (1-B)V X, pg οτε) τν... cco] (4) 

Thus 

Jie Ones) ΤΙ M NEN. 

F(x) 3x? E |n. s φίς) E ών -= 
= (1-8)v δρ ΠΕ ορ. ο(ε) | |- - B^ (5) 





From Equation (5) if is obtained that: 





2 
nen (a) 
2 
Qx 
(6) 
2 à 
(DB +E JCE) » - ló(t)*(1-8)VI,pgé(t)*pg À C(t) (0) 
Using the following relationships: 
2 D ue 1 
e E US LE DC 
a Σ (ΠΒ v(x +DB™) 
a a 
and k = Ak + 1, Equation (6b) becomes 
Ro(t)= [(1-B)Ak-B]¢it) + —*B> A c(t) (7) 
Σ ΡΒ 
a 
Let P(t) be the time behavior of the reactor power 
generation 
ο τε» οτε (8) 


f 


then P(t) € Σε b(t), and the equation describing the 


reactor power is 


€ Le 
Y + DB” 
a 


&P(t)e [(1-B)Ak-B]P(t) + a c(t) (9) 


It will be seen in the next paragraph that for certain 
kinetics problems a simplification can be obtained by 
assuming the infinite delay time in delayed neutron emis- 
sion: 


το... ο 


m 
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his expression is used as- a first estimate for the precur: 
sor concentration. 


Then 









LP(t) = [(1-B)Ak-B] P(t) + [1+4Ak] B P(o) (10) 
as: Concentration of Precursors Equation 
5 m k - E 


For a slab reactor and using one-delayed neutron group, 


Equation (11) can be expressed as 


> 


Ex, EINER) ERA) :B Σε O (x, E) (12) 


Q2 


E 
EMESUnP Equations (3) to separate variables, Equation (12) 
becomes 


Ge MC E SE E. $(t) (πιο, 


Then, using Equation (8) yields 


C(t) = - À C(t) + 


mje 


B CPC) (14) 





From Equation (14), the steady state concentration of the 


> 


precursors can be obtained as: 


org: = PO» (15) 


This result was previously used. 
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3% Energy qua cion 


Expressing H(r,t) as the energy content per unit 
volume, then the time rate of change of H must be given 
by the net energy gain per unit time and volume from the 
fission reactions and the reactor coolins: 
H(r,t) = ec [T(r,t) sur Ce) D T Om) (16) 


and 
- EET I ON o (17) 


For a slab reactor, the energy balance becomes 


oc, — 6(x,t) = el, Mao -P(x,t) (18) 
Pet OG tarot eo) (a) (19) 
Cx, tee) Cl MEER (b) 
Uc) ο Ce) F) (c) 
Then 
o O; (20) 


Let Ce = PC, and using Equation (8), Equation (20) becomes 





c 6(t) = P(t) - A(t) (23) 


Three cases are normally encountered for the power removal: 
(a) Adiabatic, @ t) = 0 
(b) Constant power removal, (P(t) = A 


Ge). Newton's ilaw of cooling, (? (t) = h$(t) 
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4. Feedback Models 
a. Temperature Feedback Model #1 
The suPttplieation factor 1S temperature and 


time dependent 


το τε (22) 
Bepandine it in powers of [T(t) - T,» yields 
eo n 
k(n) = a a, te) - ΠΡ 
n=0 
thus = 
k(T) =k [1 + >= E απ ην (28) 
n=0 
which becomes, upon dropping ali terms beyond n = 1, 
k = k [1 + ΠΤ = (1 + Ak ) (1 η (24) 
ΟΥ 
Ak(t) = Ak. eva (255 


where Ak, is the change in k applied to the reactor at 
t - 0 and Y represents the temperature coefficient of 
meactivity, which could be either (+) or (-). 
b. Feedback Model #2 

One of the major considerations to be made when 
accounting for the temperature dependence of the multipli- 
cation is that described by the nuclear Doppler effect. 

Thompson and Beckerley (23) expresses that 


dependence as dk _ 3010 η 
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here T is in ~K. 


From (26), Ak can be expressed explicitly as 


a 
. D 200m O 
Ak = Ak + qu = ) T € Το) ] (27) 
being n = 1 a special case, for which a logarithmic expres- 


Sion LS found. 


Typical váluss of nare 1/25 md 34220 


6. DISTRIBUTED-PARAMETER REACTOR SYSTEM 
E Diffusion Equation 
The diffusion equation is used in this case without 
resorting to any approximation, then 
ο ο = 1a 020) _ 
BUS (o, t) 2,0(r,t) e Est) 


-(1-B)VE ¿Pg $(r,t) - Pg 3^ À, C, Cr, t) (28) 
i 


Eor a slab reactor and using one-group delayed neutron, the 
diffusion equation takes the following form: 


9° 1 9 
Do B A QUE E 22 (x , t) = 


-(1-8)v Σε ps OO ο) (29) 


The cross sections should be written in a proper way as 
time-dependent but these variations are assumed negligible 
in comparison to the rapid transients in d, C and T. The 
equation for a non-trivial, steady state solution $o (2) > 
c, O90 dcr. 

MES 


Us 2 2.0,(x) = -(1-B)v2Z, pg 0 (x)-pg À C, (x) (30) 
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New subtraetine Equation (30) ἔτοπ Πηπαείου ο τοις 


2 
DFO) 2 2, V(x,t) = Za, t) = 
-(1-β)ν Σε ρε p(x,t) - pg AP (x,t) (31) 
where V(x,t) = (x,t) - CS (a) (32) 
= AS (b) 


This change of variables moves the equilibrium state of the 


system from (> o> Ην to MORO), Figure 3. 


ν Σε Pg 
Define k = «SE (33) 
α 
am I 
: 2. 2 
and using α - lt*L Ἑ (34) 
the diffusion equation becomes 
2 
E ee NOSE) E. 
ax? S = 
l 5 t 
+ pg A ÜÉ Guo) 21939050 (35) 


Dimensionless space and time expressions are introduced at 
HIS point: 


pos μου (a) (36) 


T = vit | (b) 





k 9 T 
+ [(1-B)ko-1]v(y,t)+ S20 AB Cy, t) 2» 2920522 | (37) 
f 
The boundary conditions are: (a) Ņ(W,t) = 0 (38) 
(b)W(-W,T) = 0 
x 
where W = Ds expresses the dimensionless half-width of the 
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Steady state 


(6,5 C.) 
V 
Steady state 
(0,0) 
Figure 3. Translation of equilibrium states. 
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cab reactor including the extrapolation distance. -Ire 
Martial conditions ame: (a) b(y,o) = 0 (39) 
(b) 2 Cy,0) = 0 
ee Concentration οἳ pzeeursors equa tmon 


dC. (r,t) 


3t = -A,C,(r,t) T vB, 20 (r,t) (40) 
EL lab reactor and using one-group delayed neutron, 


Egustion (40) ig written in the following way: 


ουν εἰ 


3c Zu Netzer vB Σε ġ$(x,t) (41) 


The same assumption used in the diffusion equation, with 
respect to the cross sections is used here. A non trivial, 
Steady state solution 9,0%); ος (κ) is assumed: 


0 = - AC Cx) + VB Ls $ (x) (42) 
EubE23cting Equation (42) from Equation (41) yields 


κ = A (x,t) + vB E, Y(x,t) (43) 


Using Equations (36), Equation (43) becomes 


vi. 2 A a) E Nee (yr et) + vB 2, Víy,T) (44) 


Mie initial condition for the equation is 


G (y,0) = 0 (45) 


Introducing the multiplication factor in Equation (44) yields 








τ) SaL k 
E πο “ο. I 0 y) 





ot 


a OT (46) 
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Figure 5. Dimensionless 
slab reactor. 
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Ss Energy Equation 


The energy balance condition. τὸ normally expressed 


as 
Ener roduced - Ener removed = Energy stored in 
ze ü Ba the system 
The energy produced is due to fission: € » $(r,t) 


The energy removed is due to the coolant: (7 (r,t) 
Three cases are normally encountered for the energy removed 
from a system: 

(a) Adiabatic case, (P (r,t) = 0 

(b) Constant power removal, for which the steady 
state requirements are normally used, (? (r,o)- ΕΣ εφίτ, ο) 

ΟΥ "UE la (i) 
(c) Newton's law of cooling, P (r,t) =h 9 (r,t) 


It is noticed that these cases are for stationary-fuel 


meaetors. The energy equation can be written as: 
AE d 2 
E. u = € Σε $(r,t) - P {ΠΕ} (47) 
Let Ce - De: be the heat capacity of the reactor system 
f 


expressed as energy per unit volume per unit temperature. 
The subscripts f and C stand for fuel and coolant respec- 
tively. Using the same assumption as used previously with 


respect to the cross section, and expressing Equation (47) 


for a slab reactor, yields 










oT, (x,t) 


sr EI NE) (48) 


35 





2e Adãabatic Case 


oT. (x,t) 


ος το . ue Dex a) (49) 


This case wili not be considered in this work because it 
represents accident conditions., 
Da Constant Power Removal 


oT, (x,t) 


οἱ --ἷς-----εΣ. φίκ,ε) - Cn (50) 


Assuming a steady-state, nontrivial solution 9,0%) it 
becomes 


0 - e E, ϕφίκ,ε) - Ce. (51) 


Subtracting Equation (51) from Equation (50) yields 


C ion = € E, P(x,t) (52) 
where Olot) = Tex, t) - Tso (X) (53) 


Introducing Equations (36), Equation (52) becomes 


96 ο 
Covi. a Cy) € Le. v(y,T) (354) 


Zeh initial condition: 
8(y,o) = 0 (55) 
c. Newton's Law of Cooling 


oT. (x,t) 
Ce SA = € Σ Oz) = h(T,- T) (56) 


Assuming a steady-state, nontrivial solution, d(x) and 


Τεο(κ), Equation (56) becomes 


en SUD = € E, v(x,t) - μη (57) 
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Note that T, and h remain constants, This 15 explained by 
means of the one effective coolant temperature assumption. 
Introducing the dimensionless expressions for space and time, 


Equation (62) becomes: 


JEN TUE. . 
ος V Y T € Σε WCy,T) h6(y,1) (58) 
ΠΠ initial condition: "Uly,C0) "0 (59) 


i‘ Feedback Models 
a. Temperature Feedback Model #1 
Through the temperature, the multiplication 
Motor is now dependent on both space and time. 
k= k(T) = k[T(x,t)})] (60) 


Ending it in powers of [T(xXjt) - T (x)] yields: 


k(T) u : τν [T(x,t) - T, G0]? 


Thus 


co A 
κ, ο ο ο IE) (61) 
nd 


which upon dropping all terms beyond n=l, reduces to: 


k= k (1 + y(T - T1,)] = k (1 + 1θ) (62) 
where y = 9k/9T = A, 
and O(xi t) = T(x To= T£) 


Mien, Equation (62) becomes 


| k=1+ Ak + yo | (63) 
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having neglected the term yAk_ 0. Ihe multiplication factor 
can now be written as: 

PST) 1 + Ak, ο δίνετε, (64) 
Ene temperature coestLelient 0f reactıvity, which could 
be either (+) or (-). 

b. Feedback Model #2 

The temperature dependence of the multiplication 

factor can be expressed in a general form, adding the con- 
tributions due to the Doppler effect, due to structural ex- 
puansrons and due to other effects, as: 

dk ln T m 

Eo M LONE Yg G2 Ta (65) 
where D and E stands for doppler and expansion respectively. 
k can be expressed explicitly, using the initial condition, 


k= k at T= T thus 
ο o 





y a a 
o D EU nu Md un 
k = 1+Ak + O In) T, GE } + 
Y T 
E IA m Ον 11 . 
ο ο COPIE UE ot (66) 
when n = 1, the third term on the right hand side becomes: 
Se In = (66a) 
Dre T 
ο 
It is noted that a, = a) = 300°K and T is usually expressed 


mm ^K. 
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Ill. STABILITY ANALYSIS 


A. POINT=KINETICS MODEL 

The reactor under consideration for this analysis has 
the following characteristics: 

fa) thermal reactor: 

(b) Homogeneous, bare, slab reactor. 

(c) One-group delayed neutron. 

(d) Constant power removal. 

(e) Stationary - full reactor. 
and normal operating conditions (no accidents) are assumed. 


The governing equations for this system are: 


τ, S ESSET (Cesare ice) (10) 


Ce (21) 


using the following feedback models: 
(a) Ak = Ak, + yo (25) 

D „300, n 300 
S X (---- 


l-n [TK O T 
O 





(b) Ak = Ak + Ja] OD 


where ΑΚ, is a positive step insertion of reactivity by 
external means, such as control rod motion. The initial 
conditions for the system are: 

(a) P(o) = (ο) 

(b) (o) 

Ce) aCe} 


0 


0, T(o) = T 

ο 
For the constant power removal, (P (t) becomes P(o) and 
Equation (21) can be written as 


CN ON (214) 
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Assuming P(o) = 1 and introducing a dimensionless power, 


Eustions (10) and (21a) can be expressed as 


τρίο) - τρ ου E Clee 


cO (t) 


P(t) 1 (21b) 
if, No Delayed Neutrons, No Reactivity Input 
The point kinetics model when 8 = O and Ak = 0 is 
reduced to 
ce =) yor (67) 
c0 =P- 1 (68) 
After some algebraic manipulations, Equation (67) can be 


expressed as 


I^ 


In P = = (69) 


a, 


t 


Mass multiplying" Equations (68) and (69) yields 


: Y : 
d C 
P -= at In P - i 66 = OQ (70) 
Siren can be written as 
Y 
d JN c A 
dc P - In P - 7 7 68°) = 0 (71) 


e” expression inside the parentheses can be called V, then 


Il 
H 
| 
H 
= 
H 
| 
ο 
cD 


V (12) 


and therefore e 0. 


mer V to be a Liapunov function, it has to fulfill the fol- 


lowing conditions: 


(a) V(P. 9) > O0 for Τοπ Ὁ 


(bX Vlo) 0 (73) 


α ΝΕ ο) O for P>0 and 9-0 


IA 
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IV of Equation (72) is modified as 


> 


De ser (74) 
2”. - 
NES scen that Equation (71) is not altered, and for a 


negative reactivity input Equation (74) becomes 


| y | 
C 2 
7 0 (75) 





VS (Re ~ In Par) + 


Noles 


Me expression (P - In - P-1) > 0 for all values of P > O. 
ΙΙ satisfies condition (73a) for all values of P and 
0 > 0, when y is negative; for a positive value of y, con- 


dition (73a) is satisfied whenever 


(P - In - P-1) > 23 ο 

Condition (73b) is clearly satisfied for the steady state 
Bu tion P = 1 and 9 = 0. Condition (130) gives V = 0 for 
both cases of positive and negative y, which means that 
asymptotic stability has been ruled out. A similar result 
was found by Ergen and Weinberg (24) using the Hamiltonian 
approach. Due to its simplicity, this case has been added 
to this work to acquaint the reader with a general view of 
mae necessary steps in the Liapunov's Method formulation. 
This case does not represent any meaningful practical 
situation. l 

Typical data for a thermal reactor is obtained from 
Solomon and Kastenberg (26), as shown in Table 1. From 


Muss data, £ is found to be 1.235 x TOSS An average 


heat capacity for liquid light water is used, 1 ale Gr 
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TABLE 1 


2 0.49780 cm 
0201 cm 


229] 10? o 


0.0064 


0.409718 cm 
1e E 


0.080 sec 
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The case of negative Y has been found to be stable 
Bor all values of P and O, but the stability is only 
neutral (V = 0). 

The case of positive y is represented in Figure 6, 
in which the system is seen to be unstable for most of the 
operating ranges of P and 0, 

AOne Group Delayed Neutron, No Reactivity lnpurt 


In this case B is included in the analysis and 


Ak = 0, the governing equations become 
ee πη | (76) 
có = P- 1 (77) 
Ak - YO (78) 


EUDSETrtution of Equation (78) in Equation (77) yields 
Peel) yO Db CUtyv9)B (79) 
cross multiplying" Equation (79) and Equation (77) gives, 


after some algebraic manipulations 


ο 
d - ΟΥ (1. ΤΕ ve) 
SO o + oo (tyo) (80) 
Thus 
B a 
<[P-InP- = £(1-8) 70+ = 9] = 5 (1+y0) u (81) 
Then 
B 

me -τη»1).- i A (1-8)0% + 9 (82) 
v= 5 (14y0) SG (83) 


It can be seen that the inclusion of the delayed neutrons 


complicates the Liapunov function, obtained by the same 
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ο. ο 


Figure 6. Stability domains for the case of no delayed 
Beutrons, no Ak input, positive temperature coefficient 
or reactivity. 
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steps as for the previous case, but this case approaches a 
more real situation. 


Two cases are analyzed: 


av Negative reactivity: y = = NA 
1 c | y | 2 Pe 
eer - In P-1) + 2 πὴ C15 + Tm" (84) 
Mu... (P-1) 
NS Cuomo (85) 


NES seen from Equation (84) that V > O for all positive 

values of P and 9, Normally (1-|y]0) > 0, then for values 
JEO P < 1, ν is less or equal than zero. Clearly V(1.0) 
= 0, then V only vanishes at the equilibrium state. It is 


concluded that asymptotic stability is obtained whenever P 


Πίπας than 1.0, and for all 9 > 0. Since Δε, = 0, P can- 


not be larger than 1.0, Figure 7 shows the domain of 
srabıility,. 

B. Positive reactivity: y = Iv | 
Equations (82) and (83) represent this case. The condition 


of asymptotic stability is obtained whenever 


B 


ς 


(P - In P-1) + rm θ 2 >, T (QN OS 


end P < 1. 
This case presents a very narrow region of stability, as 
seen in Figure 8. 
5.  One-Group Delayed Neutron, Step Reactivity Input 
The two previous cases have only theoretical inter- 
est because they do not represent a practical situation. 


The case of Ak input represents a real practical situacions 
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dV 


~= 0 
dt 
a 
i | EDEN y 
/ sasie /// / 
Sl) / 
/ 
| ARA τ 
0 9 
EN ure 7. Stability domains fer the case of delayed 


Zor Ons, no Ak input, negative temperature coeffi- 
er Of reactivity. 


P 
ly| = 0, Poe 
ν-ο 
3 IM > v=0 
2 
d 
απ =0 
S “7 
EI 
(= ON 
᾿ ins //, 
0 10 15 20 6(°C) 


Figure 8. Stability domains for the case of delayed 
neutrons, no Ak_ input, positive temperature coeffi- 
pon" mor reactivity. 
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meee start-up and shut-down, planned or unplanned motionkor 


Gomtrol rods, The governing equations for this system are: 
QP = [(1-B) (Ak + y0)-B] P + (1+Ak + y0)8 (86) 
foe Pp - 1 (87) 


Substituting Equation (87) in Equation (86) yields, after 


Simplifying terms: 


P = [Ak (1-8)-8] $ 0 - y(4x, + ΑΥΘ) + Y (1-8)? (88) 
Then 
ma . oey -geI Ê a ld Yere 
gefr-1-tar c1-8)-81 5 0] = gar +8y0)+ F(1-B) oF (89) 
Thus 
V = (P-1) + [8- Ak (1-8)] m 0 (90) 
"B E 
dem „(Ak + By0) + j (1-B)0P (91) 


For negative Y, Equation (91) remains unchanged and Equation 
(90) becomes 


ν - (Ak, -Bly| 8) - Axl 1 gor (92) 


In order to obtain a positive definite V, it is required 


that 


ς 


P + 18 -Ak_(1-8)] 7 9 > 1 (93) 


mrs condition is Satisfied whenever 


B S 
Ak. < m (94) 


for all 8 < 0, in agreement with the linear theory. To get 
a negative definite V, it is required that 


NS Clap | O 955) 


This is the additional stability condition required by the 
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Bonlinear theory. Equations (94) and (95) are specifying the 
bounds for asymptotic stability for a negative Y. It is 
noted that immediately after the step change in the multipli- 
aition factor, the neutron flux increases rapidly; this is 
referred to as the prompt jump, Lamarsh (25), and it is ex- 
pressed as 
B(1-Ak ) 
Pet) aT Pao ) (96) 
O 

Bor positive Y, it is seen from Equation (91) that V will 
always be positive definite for all values of P and O0, 
representing an unstable situation. Data from Table 1 is 
used to specify domains of stability for various Ok, inputs. 
The first bound for all these cases is provided by Equation 
(94) giving a Ακ; < 1.00625$,; in order to obtain a positive 
mn te V, for all P > l and 0 > 0. It is clearly seen 
that P has to be greater than one, the steady state power, 
after a Ακ; insertion. The second bound is obtained, set- 
ting V = 0 in Equation (92), Figures 9 and 10. After the 
EuSertion, the power will increase to an amount specified 
by Equation (96). From Figures 9 and 10, it can be seen 
that right after the Ak insertion, the system is unstable, 
and P and O increase, until the boundary V-0 ls crossed, 
then the system becomes asymptotically stable and the tra- 
mectory returns to the equilibrium position, but when V=0 
is crossed again, the system turns to be unstable, and in 
that way, the trajectory keeps oscillating along the line 


V=0. A similar situation is obtained by Meghreblian and 
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mure 9. Stability domains for the case of delayed 
Mmeutrons, step Ak = 0.156$, negative temperature 


coefficient of reactivity. 
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Ak, Ξ 0. σος 
4, πμ 10”*/°c 
Prompt Jump 


m nn nun a nn 
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Heure 10. Stability domains Iorlthe case ot delayed 
neutrons, step Ακ; = 0.3125$, negative temperature 
@eefficient of reactivity. i 
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Holmes (22) when solving the system of differential equa- 
tions. Their results showed that this is the case of damped 
oscillation, the presence of the delayed neutrons providing 
ο ἁσπρίηρ force, thus, following an insertion of reacti- 
vity, the system will reach a new steady state power. 
4. One-Group Delayed Neutron, Ramp Reactivity Input 

The governing equations for the system are Equations 
(86) and (87), with the difference that now Ak, is a func- 
tion of t expressed as 


Ak (t) = at (97) 


where a represents the rate of insertion of reactivity. 
Combining Equations (86) and (87) yields, after some alge- 
Ὁ "πο manipulations 


TES ui: a ορ το 6 (98) 


Then 
Y β 
d T ο ο E s 
αι ο ο ορ μωβ μα ο ua Su 
= (1-8) SP + D NE (99) 








B 
er) a-g) YE g^ uS 
V = P 1 > (1-8) 2 0 Pa 9 (100) 
ee yo Bat 
V (1-8) DEUM c 7 (101) 
For a negative Y, Equations (100) and (101) become 
2 
τι. τ πε ΕΞ: (102) 
£ 2 2 
y = (1-8) 2È p - YM y Bat 
V (1-8) o DENT (103) 





mae following condition is obtained for asymptotic 


stability: 


O 


» ue inu AEn 
Ak (t) e: mue CIE) PFE 


(104) 


Tt can be seen that V is positive definite for all P > 1 
med > 0. Equation (104) is necessary to obtain a negative 
definite V to ensure asymptotic stability. For a positive 
Y, ple »eystem is inherently unstable. Equation (101) gives 
a positive definite V του all values of ( Wand) 0 πο ο. 

A typical ramp insertion of reactivity is shown in Figure A. 
During the reactivity insertion, the system is unstable, 
then V = 0 is obtained after the ramp insertion has ended 
(10 sec) and the power and temperature have risen to an 
appropriate level. Setting V = ο in phat tom  ClO3). the 


following expression is obtained 


as BA. 
x ae US (103a) 
ο 
ie end of the Ακ, insertion (t = 10 sec), this expres- 


Sion implies that P depends linearly on 0, Figure 11. 


E Üncecnoup Delayed Neutron, CeneraliReactivity Input 


The equations for this system are: 








&P = [(1-B)Ak-B]P + (1+Ak)8 (106) 
πρ η (107) 
Yp 300.n 300.n 
A O eee τ. ] (108) 
Let 
N as 
SE ITE N (109) 
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Figure A 
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Ak _=0.1565(10 sec) 
yl 10 *yec 














va 


ΑΘ; 
AS Cable (V < 0) 





0 10 100 1000 0076) 


Mesure 11. Stability domain at the end of a ramp 
Eusertion of reactivity. 
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Then Equation (106) becomes 


1» - [(20-8) (Ak +y,*R(T)-8]P + [1+y *R(T)IB (110) 


Combining Equations (110) and (107) yields, after some 


algebraic manipulations 





LP - (1-B)Ak c0+Bc0 = [(1-8)P+8]y,R+Ak, (111) 
Then 
B-(1-B)Ak 
d ο z 
nc DE CENE pex ee 69 m 
Yp’ R Ak, 
= [(1-8)P +8] i dE (112) 


Therefore 


V 


P-1 + [B - (1-8) Ak ] e (113) 


and 
Yp'R Ak | 
ele) Bae | mm ^ (114) 


<j 
N 





Equation (113) is independent of y; therefore, V is positive 
definite whenever conditions (93) and (94) are satisfied. 


The similarity with the step insertion of reactivity can be 





seen. Equation (114) becomes, for negative y 
Ak, lern) 
A A — (115) 


Equation (115) is negative definite whenever 





Ae O (116) 
for all n # 1. The case n = 1 is the logarithmic case and 
Equation (116) becomes 

Ak, < 300[8+(1-8)P] |y,| In S— (117) 

ο 
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The case of positive y leads to an unstable situation due 


uNXUhe fact that Equation (114) is positive definite “for 


EJ n, 


P and T. Data for this case are found in Thompson 


and Beckerley (22) and are typical of a fast reactor. The 


following specific case is analyzed: 


(a) 
(b) 
Ge) 
(d) 
(e) 
(£) 


Oxide reactor with volume ratio v0, to ro, = 7 


Sodium density = 50% 


ee 1077 /*k 


n =.0,96 
B 9050053 
nl u0- 15065 
O 


The domain of stability is shown in Figure 12, where the 


V = 0 is obtained from Equation (115). 


B. DISTRIBUTED-PARAMETER REACTOR SYSTEM 


The reactor under consideration for this analysis has 


the following characteristics: 


(a) 
(b) 
NOD 
(d) 
(e) 


fast reactor 

homogeneous, bare, slab reactor 
one-group delayed neutron 
Newton's law of cooling 


stationary-fuel reactor 


~ 


and normal operating conditions (no accidents) are assumed. 
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The governing equations are: 


2 
AMO «[(-8)ka-1]9 (y) * $$ AC, - 


ΟΥ vo 
m Il (ap 
Bo 2K 
. 3€ 
vi. LL. -A-f (y, 1) ji: E s UCy, 1) (46) 
ENT. LD E Lo Uyt) nO.) (57) 
f a OT f y» y» 


using the folloving feedback models: 


E  k(y,1) -1 + ΕΙ μιν. τ) (63) 


Y a a 
E ο [ης ΩΣ. π ΣΙ 4 
ο τα m T ] + 


T 
Om 
E cL lost μα 35) (66) 


where Ακ. is an external positive reactivity insertion. The 
boundary conditions for the system are: 

(a) v(+ W, 1) = 0 
The initial conditions are: 


(a) vY(y,o) = O 


(b) (y.o) = 0 
(c) θ(ν,ο) = 0 
The case n » 1 will give the logarithmic term in the doppler 


expression. 


L. .-One=Group Delay ede Newt non yo ee prc Oey aa amano 
To formulate the concept of stability a distance be- 


tween the equilibrium position (oo; Cv. T) and any perturbed 


O 
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Ask cion (4, €, T) is defined. The variation of this norm, 
Which is induced by the inner product of the Hilbert space 
in which the solutions of the system are defined, with time 


DES provide information concerning the perturbed position. 


Let W 
2 202 25-02 242 
ο sc, 0 1,1 10 = á uU D (118) 
: € 
-W 
where , ; 
fuse 
2 f 
Ay = > (119) 
ED 
2 a 
A, room (120) 
2. ο 
2 μα. 
Az = Ce V 2 (1299 
and v= > = or 
= C= ο, 
9 --— T 
O 
The Liapunov function is then defined as 
7 E 2 ος T 
vQ,6, 9 = |lal|” =<a,d> (122) 
Dich by definition is positive definite. It is noticed 


A, and A, makes 


that the introduction of the constants Ajo 9 3 


V represent the total rate of energy increase per unit 
volume of the system, which is defined a the distance be- 
tween the equilibrium and the perturbed states. 

V vanishes only at the equilibrium position ($55 Co 


Ti therefore, the requirements of a negative definite V 


will provide the domain of asymptotic stability. 


E 





2 ( 
ay [oe Σ μην + Ea e 3€ “5 Es 29 C dy 
p? a dy 


n ΟΥ T 
-W (123) 


w ti Cuting Equations (37), (46) and (57) in Equation (123) 


yields 
dV w, Κλα 
dc" d ΠΣ + [(1-ß)ka- 1]v^ + νε έν + 
2 2 
vh ka YAA ,2 20.8, κα 2C,v2 h 2 
T : 2 ae Sy 2 3 2 us 2 5 s 2» 
Bv Σε B^v Σε νεΣε ε Le 
(124) 
Integration by parts of the first term and using the 
boundary condition Y(+ W, T) = 0 yields 
W W | 
9° dp, 2 
] Y m dy = -| dy | (125) 
dy? PEL. 
-W -W 


Knops and Wilkes (27), in their work, provide the inequali- 
ties useful for this kind of analysis. The so-called 
Beueenvalue inequality" plays an important role in this 
specific problem. For a system defined by the eigenvalue 
problem 

Vu + Au = 0 
in domain with u(a,t) = ul-a,t) = 0, on the boundary, the 


following inequality can be stated: 


λ / μ᾽ γς / E dy (126) 
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En inequality that can be proved using the calculus of 

ta tions. (See Appendix A). In applying this inequality 
to this work, it is assumed that the perturbed reactor has 
the fundamental eigenvalue not appreciably different from 


imac of the unperturbed reactor, then 





\= —— = BYL (127) 


Thus, Equation (126) becomes 


W W 
jas dy > BYL J y? ay (128) 
~W -W 
Substitution of this result in Equation (124) yields 
dV p 27927 2 
= - ε2Σ BOLO» - [(i-8)ka-1]v^ - 
dt — f 
ZW 2 
ιτ τη và À 
κλα 2 dd a 2 
ey —$—— pe A _ E 
Ve E: βν2Σ 2 é β2ν2Σ 2 6 
f í 
2 
2C,v2 ka 2C vi h 
E gy 4 ——— q? dy (129) 
vez. ο... 


61 





ο. ον introduction of the feedback term, the following 


expression is finally obtained: 


W 
αν WEE 2 Br , 2 
me rel, a f E (1-8) (1+4k 3] y 
i -W λ v 
* (1-80y8U" - (1I+Ak ) Co— +t TD Y ib 
T Bv Σε 
λ T vb A 9 
ντ, a 
2 
2C.v2,h 2Cevb, 
+ E . 9 ps a) ER 
€ Le α EL, y 
en. 
He 03 αν 
2» ν 
Pet a (1-8) (1+4k_) 
ULLA 
a, = = 
Z ο ο 
ne meee 
3 20 2 
€ Le a 
v : 
À a 
a NEU 
f ENS 
f 
2 
uM 2C.v2 
> 2 
EL, 
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(130) 


(131) 


(132) 


(133) 


(134) 


(135) 





Bien Equation (130) can be expressed as 


W 
dV - 2 2 2 - 2 2 
Ee ε Le α | Jis W t a, P + a, $^ 
-W 


(1-8) y0p? - (1+Ak ) a, a + a, YOY - 


-- 


(IHK) a, OU + a, γθ2ψ] dy (136) 


To obtain asymptotic stabi'trty, dt &senecescsary to prove 
that the expression inside the curly brackets is positive 
definite. A strategy successfully employed by Buis and 
Vogt (28) will be used here. Equation (136) can be written 


as follows: 


W 
dV E 2 2 2 ον 2 : 
RN e. a Cf mn + a, É + az 0%)dy] 
-W 
W 
J itak )a,Ev=a,YOEY-(1-8)10Y?+(14Dk )a¿Oy=a¿Y0%p]dy 
O A A 
W 
2 2 2 
Ja. U^ + a, © + a, ο (137) 
-W 
Equation (137) can be expressed concisely as 
dV 
dt Se mr P Q (138) 
where R - eI, a is essentially a positive constant, 


W 
a 2 2 2 
P = / (a, o^ sb aot * a. 6“) is greater than zero 


whenever aj 82 and a, are greater than zero. It is already 


known that 82 and a, are positive constants. 


Thus a, = 1 - (1-8) (1+Ak ) has to be greater than 


Zero in order to have P ^» 0, 
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απ this condition, the first bound. for ¡stability 15 


obtained 


(139 





Ia e reement with the*"linear theory. To prove that Q, 
second expression in brackets in Equation (137), is positive 
definite, the Buniakowsky-Schwarz inequality will repeatedly 


be used 
τ ο ορ) ο. ο 
f sg ay « C f A (140) 
-8 -8 -a 
Due to the positive Ak, inse®#t#6n, the syste@'s state 
variables ψ, é, O are positive for any t > 0; therefore, 
the absolute value in the inequality is not necessary. It 


is also known that 


fa? m d pa]? (141) 
-W 
ER 2 2 a2 
y Aj; ay s [3] (142) 
ο. ED 
ds of ay < [al (143) 


State variables in engineering have the properties of func- 
tions in the Hilbert space, and the norm in this space is 


specifically defined as 
a 


Cf v^ax) ^? = hid (144) 


-a 
Due to the fact that the state variables in the present 


work, v, É and 0, are continuous in space and time domains 


together with their derivatives to arbitrary kth order, it 
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is expected that they also have the properties of functions 
in the Banach space. The norm in the Banach space CL? can 


be defined as: 


a 
( {ν᾽ ne - Jul, for mo l (144a) 


m 
Then, repeated applications of these inequalities yields 


(1+Ak ) a 


ο jm 2 ek (1-28) Mi m emm 
τω ie Dal πμ + ο 4} 
mm | να. 7 AJA, N 
es ME Uo A o9 
a a a 
min (= : -—, =, fall’ (145) 
ie ek ER 


where the minimum is taken in order to assure the inequality. 


Therefore, for Q > 0, it is necessary that 


(1+Ak_) a, as | El 82 a, 
(+ 2) - min (4,, 4,, 4) 
A A A 2 2 2 
E 1 Z 3 1 A, A, 
lal] > : : (146) 
2WY ( 4 y (1-8) O ) 
A ALA Α.Α 2 
1 2 3 173 A, 


This is the additional stability condition prescribed by the 
nonlinear theory. The inequality (146) gives the second 
bound for asymptotic stability of the system, which repre- 
memes a spherical surface in the first octant, i.e. Y, É 

and 06 > 0. Domains of stability are represented in Figure 
13. The solutions for É (y,1) and O(y,1) can be obtained 


irom Equations (46) and (57), as follows: 


T | 
É (y, Tt) = η τ ας (147) 


ο 
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W 


cay" J goa 


-W 


Stable 


(Ν 


«C Trajectory 


1., 0 -{0) Lote | 
x m > c c cay? f ψ αν) 
Unstable Es 
W 
a f 924,5 1/? 
-W 


BXoure 13. Schematic of the spherical surface determin- 
ing stability domains for the distributed parameter 
reactor system after a Ak. insertion. 
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where A (τ-τ') 








=r! = 
O min vz ¿8 e (148) 
and 
E 
B(y,T) = [nem UT Pd (149) 
Ó 
where ST) 
Ed; Cuán : 
H4(t-1') e (150) 
Ον δ. 


Application of the mean value theorem, in the time domain, 





yields T 
L 
Ἔ(ν,τ) - νο, Ὁ / "πιο (T51) 
ο 
where λ 
βννΣ ϱΣ, mos τσ 
Τη (τ) AA (1 - e a ) (155.2. 


aad «T x T., 


Also: 
E 
στ - νο,» / E O dr Vp T. (x) (153) 
where 
_ Er z E E / 
το (1) πες (1 - ο a ) (154) 


and O « T< T. 
EUbSstituting Equations (151) and (153) in Equation (118) 


results in 
W 


W 
pa = J ODay +f A, 92 (y, D) T, ^ CO dy + 


- -W 
W 


+ favo DES (155) 


- Ww 
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Applying the mean value theorem, in the space domain, yields: 


I 


D = a, + 2,9 GD 1,500) 


* 24A, V^ (y, 1). T,^ (1) (156) 
or 
BU = y(y,7) ο. A D GE) 


where O « t « r-and -W « y ob... 
Then the mean value flux can be expressed explicitly as 


follows: 
E" |a | 
Mam = (158) 


Y RII a) ο 


2 2 
1 2 Τρ” (τ) 


3 
Typical values for the fast reactor nuclear parameters are 


listed in Table 2, after Solomon and Kastenberg (26). Using 


these parameters, the constants for the system are found to 


be: 
a, = So UM sS πο am sea 
a, = δ 155.55 D” En ON 
a, = ο ο ο 10° en-sec 
a, = O xc > ο n 
a, e η. O μάς. 
D 2.780 x or T 
= 1.920 x m a E EA -sec 
W= 3.0 
and 
T,(t)= 70.05 (l-e N (159) 


12,, ,-0.00004481 


SPEI 2 ΘΕΣ ET (I. ) (160) 


68 





Dom Equation (139), the first bound is given by: 
Ακ, < 22.003318 

Brom Equation (146), the norm, second bound, is obtained 
to be: 

|i] » 2 x10 
for a ΔΕ, s 0.135655 The mean value flux from Equation, (loc 
sehen in Figure 14. „The physical meaning of this mean 
value flux U ο) is depicted in Figure 15. The mean value 
flux provides an order of magnitude information of the in- 
crease in neutron population at the time the surface of 
Stability is reached. From Figure 14, it is seen that a 
fast rising flux will quickly reach the surface, whereas a 
more slowly-rising flux will take longer to reach the sta- 


bility domain. 


ELE Ome-Group Delayed Neutron, Space Dependent 
aep Reactivity Input 


This system is a special case of the previously 
analyzed system. It is desired to find an answer to the 
following question: Given a positive reactivity insertion 
Ακ. is it safer to insert this "hc uniformly across 
the reactor, Figure 16a, or to insert less reactivity in the 
central region, Figure 16b, provided the total reactivity 


insertion is the same? 


Equation (129) can be rewritten as: 





W vÈ À 
PAR πμ 
-W B^v Le Q 
2C¿vE_h a DES 2C,v2 = 
emer O° (λε T ὁ ν - —À———— ow dy (161) 
εὖΣεα f m ve,” | 
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Figure 14. Mean value flux vs. dimensionless time plot. 
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Figure 15. Mean value flux schematic in the slab 
reactor. 
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Figure 16. Space-dependent Ok insertion: 
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Mation (161) can be rewritten as: 


-2W/3 2W/3 W 
ST < -era I f {-}dy + f i-}ay+f -} αν (162) 
=W -2W/3 2W/3 


Three steps Ak, are introduced in the following way: 


k = 1 + Ak, -~ YO τος MOS yE e zt (163) 
i > E ο d 2W 
k = 1 + Ak, yo for 3 u 3 a) 
29 
k = 1 + Ak, - YO τος 3- < y <W (165) 


EbStituting Equations (163), (164) and (165) in Equation 
(163) yields, after some algebraic manipulations and intro- 


X : 
duction of the constants 81» 85» Az» 8, and as 


-2W/3 
dV _ e?y 2 12 2 2 u 2 
ga TEL A [a, yor ant + a,0°+ (1-8) yep 
-W 
~(14Ak,)a, G+ a, yOE j- (1-Ak,)a,0y * a,.y0^y]dy^ 
ae ayot ES 5 


2W/3 
+ of han yt aj € ^ a 40^ (1-8) Y0U^- (14Ak, ) a, C V + 
-2W/3 


" a, YO (V - (1+4k,)a¿0y + a¿y0*p] dy + 


af DE ajQ “+ a, 0%+ (1-8) yOW?- (144k, )a,€ V 2 


2/3 
4 a YOG y -(1+4k,)as 0y + az y0*pldy (166) 
where 
m = 1 - (1-8) (1+4k,) (167) 
ΟΡ» (168) 


zs 





Following the same steps as in the previous analysis, the 


following conditions for asymptotic stability are obtained: 


B 


B 
Ak, yd πρ (169b) 


conditions which agree with the linear theory, and 


(1+Ak,) a 
































ο ος 1 2 3 
xc d Mi oum ; : ) 
i A i AG "€ ks ΝΕ 
[417 Ξ : (170) 
lc t a n 23) 
(i 5299 ES A. 
(1+Ak,) a a a a 
ee | A a) 
7 1 2 3 A,” Awe A. 
lal, > : - (171) 
| a Let y a i 25) 
jo πα. 


These are the additional condition of stability required by 
the nonlinear theory. In order for the total reactivity 
insertion to remain the same for both strategies, one must 


have 


-W -2W/3 2W/3 

1 X EE 

kJ Ak dy = = J ia i Ak, dy (172) 
W E AD IG 


Pom Figure 17, it is seen that in Region 2 in which Ak, < 
ΔΚ, all, > | ἆ |], so that the region of stability is attained 
sooner than in the case of uniform Ak, insertion; but |4||17 
lal, so in Region l of the reactor, where Ak, > ΔΚ» the 
region of stability is attained later than for the case of 


uniform Ak insertion. However, the flux in the central 
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Figure 17. Spherical surfaces for the case of space- 
dependent Ak, insertion. 


τω 





zen of the reactor plays a more important role in the 
transient behavior of the reactor; consequently, the space 
dependent reactivity insertion, with Ak, < AK is safer 
than the uniform Ak. insertion. 

3. One-Group Delayed Neutron, General Reactivity Input 

For this case, Equations (37), (46), (57) andi (Ce) 

ape the governing equations. “The “feedback -coeffiedemes, y, 
Yo and Yg» are assumed to be negatives; Yp and Y, are the 
Doppler and Expansion coefficients of reactivity, respective- 
Ihe same Liapunov function, Equation (122), is used 
and a similar development previously discussed leads.to the 
IlENNO ine equation, in analogy do Equation (136): 


dV 


~— = 


W 
Dux emet f cuts 2,6 + ren: 
+ a, ΨΕ(ϐ) αςθψΕ(θ) - (HAK day v- (148k) aso ay} (173) 


where the Feeab#r termgis: 


iR - je |Yel?" 34m eee 
Ε(θ) = ly lo ΠΕ ΗΝ [s - το j+ ml! mec ] (174) 
Mere R and P are constants equal to 300°K. F(6) is positive 


Er all values of n > 0 and m > QO. It is noted that the 
Doppler feedback is the primary feedback mechanism in large, 
ceramic fueled fast reactors, and ans on is usually the 
dominant feedback in small metal fueled fast reactors. For 
BI, Equation (66a) is used. A different approach, more con- 
servative, is used to analyze Equation (173). The expression 


inside the curly brackets in Equation (173) has to be posi- 


tive definite in order to ensure asymptotic stability. With 
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this in mind, a comparison of the orders of magnitude among 
the terms is carried out. 


W 


Ege - o ‘| [a V t É taS - (+Ak Jay} + 


+ 8 (a0 i (1*Ak )agU] + VF CO) C(1- 8) ja f a¿0)]dy (175) 


In order to get a positive definite-furctóben inside the 


curly brackets, the following conditions have to be satis- 





fied: 
$ ς β 
(a) a, > O which yields Ak. < 1-8 
(b) ve (1+Ak |) a,V > 0 which yields 
E: E 
É (1*5k,)a, 
and 


(c) a 40 - (1+Ak_)acv > 0 which yields 
ISBN 
0 (1+Ak )a, 


Typical values for the feedback coefficients are: 


O por, Evaluating 


y= D re Ks 7x10 9/*x, Yg7 or 
these conditions for a typical fast reactor, with a ΔΕ, ~ 
UNIDOS, it is obtained that: 

(a) Ak. <2 OU Soe 

(b) PE « 3634.0 

(c) $ « 4.203 x 10% 
These three conditions have to be satisfied simultaneously. 
The domain of stability is shown in Figures 18a, 18b and 
l8c. It can be seen that the domain of stability is very 


large when the coefficients of reactivity are negative. 
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Eure 18a. Stability domains for the distributed 
Parameter reactor system after a general reactivity 
insertion. 
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Figure 18b. Stability domains for the distributed 
parameter reactor system after a general reactivity 
insertion. 
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mesure l8c. Stability domains for the distributed 
Parameter reactor system after a general reactivity 
insertion. 
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νι CONCLUSIONS 


Liapunov's Second Method has proven to be a useful 
method to find stability domains in Nuclear Reactor Control. 
The fact that a Liapunov Function is not unique makes the 
choice of such a function to depend upon the experience and 
ingenuity of the analyst. Once an appropriate choice of the 
Liapunov Function has been made practical results can be 
found. 

For the Point-Kinetics Model, a Liapunov Function vas 
obtained without imposing an a priori bound, as most of 
the conventional methods do. Domains of stability are pre- 
sented for different cases using typical data from thermal 
meactors. 

The Distributed-Parameter Reactor System was studied 
mene the concept of a norm, an a priori bound, and the 
pmarration with time of this norm, defined in a Hilbert 
space, was analyzed to obtain a domain of stability. Data 
from a typical fast reactor was used to define the spherical 
surface of stability and to estimate a mean value flux, in 
Space and time, when this surface is reached. The case of 
Space dependent Ak. insertion, with less reactivity in the 
central region, was found to be safer than the case of 
uniform Ak, input. 

The linear theory stability condition appeared through- 


out the analyses giving confidence to the obtained results, 
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ee ddition to that condition and additional conditione) 
appeared in each case, required by the nonlinear theory. 
The inclusion of the multigroup delayed neutrons is 


recommended for further improvement of this analysis. 
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APPENDIX A 


THE EIGENVALUE INEQUALITY 


Let à > O be the least eigenvalue of the system 


2 
Veen (A1) 
dy* 

with boundary conditions: (a) y(+W) = 0 


the following inequality can be stated: 


W id 
af vay Es J G dy (A2) 
a 


-W 


Proof of this inequality can be done using the calculus 


Gi variations. 


Let | 
W/ 
Je? ayhay > 0 (A3) 
-W 
d 
Applying Euler's equation: — F ,= F A4 
pplying Eu q Agee y (A4) 


Equation (Al) is obtained. 


Solving the eigenvalue problem X is obtained to be: 


A = στο = BL (A5) 
4M 
Then W W 
ΜΓ > ar yay (46) 
- WM ay 
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